Invasion-wave induced first-order phase transition in systems of active particles 
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An Enskog-like kinetic equation for self-propelled particles is solved numerically. I study a density 
instability near the transition to collective motion and find that while hydrodynamics breaks down, 
the kinetic approach leads to soliton-like supersonic waves with steep leading kinks and Knudsen 
numbers of order one. These waves show hysteresis, modify the transition threshold and lead to 
an abrupt jump of the global order parameter if the noise level is changed. Thus they provide a 
mechanism to change the second-order character of the phase transition to first order. 



Collective motion of self-propelled agents is a key fea- 
ture of active matter systems and has attracted much 
attention [l|-|3|. Systems of interest range from animal 
flocks [4] , to human crowds [5] , actin networks driven by 
molecular motors [6j , interacting robots |7j , and mixtures 
of robots and living species [8| . 

Most of our theoretical understanding of collective mo- 
tion comes from two sources: (i) computational studies 
of particle-based models and (ii) phenomenological trans- 
port equations which are usually postulated by means of 
symmetry arguments as in the seminal work by Toner 
and Tu [9|, [lfj. These authors showed that even in a 
two-dimensional system, long-range orientational order is 
possible due to the nonzero speed of t he p articles. Many 
of the computational approaches 
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the minimal Vicsek- model (VM) [2lJ. In the VM, point- 
like particles move at fixed speed and try to align locally 
with their neighbors but do not succeed completely due 
to the presence of some noise. As the noise amplitude de- 
creases, the system experiences a phase transition from 
a disordered state, in which the particles have no pref- 
ered global direction, to an ordered state, in which, on 
average, the particles move in the same direction. 

This paper is based on a third angle of attack - the 
kinetic theory of gases. Apart from a few exceptions 
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24| the kinetic approach to active particle systems 
has been less popular. However, it is very powerful and 
allowed the solution of a long standing problem, the rig- 
orous derivation of the hydrodynamic theory for the VM 
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32h34| for alternative derivations. 

revealed that 
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Direct simulations of the VM 
right at the onset of ordered motion, large density waves 
emerge. It has been intensely debated Ill4l3l l2ll 1361 ] 
whether this order-disorder transition is continuous or 
discontinuous. By now it is generally accepted that at 
high particle speeds the transition is discontinuous with 
strong finite size e ffects. T his is in puzzling contrast to 
mean-field theories 
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■ 28( which should be valid at large 
speeds but predict a continuous transition. 

In this Letter, I show that the solution of this puzzle 
lies beyond hydrodynamic theory. I find that a special 
soliton-like density wave, which can be analyzed by ki- 
netic theory but not by hydrodynamics, is able to alter 
the character of the phase transition from continuous to 



discontinuous. I calculate the global order parameter for 
collective motion and show explicitly how its hysteresis 
and its unique finite size effects are related to the prop- 
erties of the density waves. 

The reason why these waves escape hydrodynamic 
treatment is that they violate a basic principle for the 
validity of a hydrodynamic theory - the smallness of the 
Knudsen number - which requires that the average dis- 
tance particles travel between collisions is much smaller 
than the length over which hydrodynamic fields change 
considerably. This is not true for the waves emerging 
near the onset of collective motion. They are so steep 
that their Knudsen number is always of order one. 

In one of the first analytical studies of active particles, 
Bertin et al. 24J, |25| have also analyzed soliton waves. 
However, the calculated density profiles (Fig. 7 of Ref. 
J25| ) bear little resemblence with the actual profiles ob- 
tained in direct simulations of the VM, which have a very 
sharp leading edge. Gopinath et al. [30] also calculated 
waves which look different from the ones observed in sim- 
ulations. Both groups obtained waves within the hydro- 
dynamic approach and did not observe that the waves 
have any effect on the order-disorder threshold. The 
waves calculated in this Letter by means of kinetic theory 
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are qualitatively different from the ones of Refs. 
because (i) they shift the transition threshold and mod- 
ify the character of the flocking transition from second 
to first-order, and (ii) their profile semi-quantitati yely 
agrees with the ones measured in direct simulations [31| . 

A first clue about the inadequacy of hydrodynamic 
equations for the VM comes from Ref. [26j where it 
was shown that if all coefficient in these equations are 
rigorously derived from the microscopic dynamics, long 
wavelength density modulations evolve into waves with 
infinite amplitudes. Thus, contrary to Refs. 25, 30] no 
solitons could be found. The equations were derived un- 
der the assumption that higher order gradient terms are 
negligible which is not justified when steep spatial per- 
turbations of a homogeneous state grow sufficiently large. 
Usually, such perturbations are stabilized by higher order 
nonlinear terms. This is not the case here, the hydrody- 
namic equations are driven out of their range of validity. 
To discover the final fate of these waves within the hydro- 
dynamic approach one would have to explicitly sum gra- 



dient terms of all orders, which is practically impossible. 
I circumvent this obstacle by abandoning hydrodynamics 
entirely Instead, I numerically solve the space and time- 
dependent kinetic equations of the VM. Because this does 
not involve any gradient expansion, a summation to all 
orders is achieved implicitly. 

In the VM, N pointlike particles with positions Xj(i) 
and velocities v,(t) undergo a discrete-time dynamics 
with time step r. The evolution consists of two steps: 
streaming, where all positions are updated according to 
Xj(i + t) = Xi(t) + TVi(t), and collision. The magni- 
tude vq of the particle velocities is kept constant, only 
the directions Qi of the velocity vectors are modified in 
the collision step: a circle of radius R is drawn around 
a given particle and the average direction <&i of motion 
of the particles within the circle is determined accord- 
ing to $i = arctan[^"sin(0j)/^™cos(0j)]. The new 
directions follow as 6i(t + r) = $j(t) + £». Here, £j is 
a random number which is uniformly distributed in the 
interval \—T}/2,r]/2]. 

Following Ref . |26j , the time evolution of the VM can 
be described by a Markov chain for the N-particle prob- 
ability density. This equation is exact but intractable 
without simplification. The easiest way to proceed is to 
make Boltzmann's molecular chaos approximation and 
assume that the particles are uncorrelated prior to a col- 
lision, which amounts to a factorization of the N-particlc 
probability into a product of one-particle probabilities. 
Because this assumption neglects correlations and leads 
to an effective one-particle picture, it can be thought of 
as a sort of mean-field theory which looks like an Enskog 
equation, 
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The distribution function /(#, x) is proportional to the 
probability to find a particle with a given angle 8 at lo- 
cation x. Details on this derivation can be found in Refs. 
[H [M S3|- The r.h.s. of Eq. {!]) is the collision inte- 
gral and will be denoted as /[/]. It is a nonlinear func- 
tional of the distribution function with a singular kernel 
which consists of the periodically continued Dirac-delta 
function, 6(a) = J2m=-oo ^( a + 27rm). The argument of 
the exponential in Eq. (J]), Mr(x',£) = J p(y,t)dy, is 
the average number of particles in a circle of radius R 
centered around x' = x — vr where v = vo( c os(9,sm#) 
is a velocity vector. This interaction circle is not cen- 
tered around the final position x because after the refer- 
ence particle i = 1 has collided with particles 2, 3 ... n 
it will be convected to location x in the subsequent 
streaming step. The symbol denotes spatial integra- 



tions over the collision circle. The particle density p is 
given as the zeroth moment of the distribution function, 

p(x,t) = j* n f(e,x,t)de. 

Eq. (HJ) is solved by an algorithm which is related to 
the Lattice-Boltzmann method (LB) 41]. It relies on a 
set of Q microscopic velocities, e^, where every velocity is 
associated with a distribution function /j(x, i). The posi- 
tions x are discretized on a regular lattice of size L x x L y . 
See Supplemental Material at [URL] for details on this 
numerical method. In contrast to LB a very large num- 
ber of velocities, Q « 1000, is used to accurately resolve 
the order-disorder transition. Another difference is the 
nonlocal collision term /[/] which requires spatial and 
angular integrations. Naive attempts to calculate I[f] 
by a simple integration scheme lead to prohibitively slow 
performance. A much more accuate and faster way is to 
evaluate the collision integral in angular Fourier space. 
Expanding /, 
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performing a similar expansion for I[f], and evaluating 
Eq. ([1]) in this basis gives a simple set of algebraic rela- 
tions for the Fourier coefficients of the collision integral in 
terms of gt- and hk, see Eq. (S2) in the Supplemental Ma- 
terial at [URL] . Three-particle and higher order collisions 
have been neglected. While this restriction to binary in- 
teractions keeps the simulation times short and reduces 
the validity of the numerics to low densities, Mr < 1, 
it is not a principal limitation. Similar to Ref. [35], 
the algorithm can be easily extended to include genuine 
three-, and higher particle collisions. In Eq. (J5|) all an- 
gular modes with wavenumbers k larger than the cut-off 
value kc — 5 were neglected. The remaining modes are 
sufficient to describe the behavior of the order parameter 
in the vicinity of a phase transition. The global order pa- 
rameter O is defined by means of the k — 1 Fourier coeffi- 
cients, fl = (\J g\ + hf) where the brackets denote an av- 
erage over the simulation box. These coefficients are pro- 
portional to the components of the momentum density w, 
g\ ex w x , hi oc w y , which is given by the first moment of 
the distribution function, w(x, t) = J Q v(0)f(8, x, t) d0. 
Let us first consider a very small system, L x = L y = 4 
with periodic boundary conditions and average parti- 
cle density p = N/(L x L y ) = 0.00424. To initialize a 
disordered state, all angular Fourier coefficients except 
.9o = Po/27r are set to zero. In this and all following 
simulations, the average particle number in the collision 
circle, M = irR 2 po, is set to M — 0.03 and the time 
step is r = 1. Recently [26(, the critical noise rjc a s a 
function of M has been calculated. By choosing a noise 
value r\ — 0.43 slightly smaller than the threshold value 
r/c = 0.4361, the system is quenched into the ordered 
state. Since the system size is much smaller than the 
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FIG. 1. (a) Steady state order parameter Q; (b) speed of 
the invasion wave vw (circles) and average particle speed 
u = |w|/p (triangles) measured at the tip of the wave versus 
system size. The dash-dotted line shows the speed of sound 
in the disordered phase, vs = «o/v2. All speeds are plotted 
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FIG. 2. (a) The Fourier coefficients gt for a stationary inva- 
sion wave travelling into the positive x-direction as a function 
of position, (b) The rescaled density ps — pLZ 2 10 6 versus 
rescaled position Xs = xL 2 x lQ~ b for different system sizes. 
Parameters are the same as in Fig. [T] 

critical length L = 2ir/k (see Fig. 2 of [26|) above 
which the homogeneous ordered state is linearly unsta- 
ble, the system is expected to stay homogeneous. The 
numerical solution of Eq. ([1]) agrees with these predic- 
tions: the zero momentum disordered state evolves into 
a stable ordered state. As shown in Fig. QTa), the order 
parameter stays constant at increasing system size until 
a critical size of L* ss 270.5 is reached. When the system 
size is adjusted by just one lattice unit from L — 270 to 
271, Q makes an abrupt jump, almost tripling its value. 
A closer look reveals that while the steady state solu- 
tion is homogeneous at small L, it is not homogeneous 
above L*. Instead, a single-peaked density wave is go- 
ing through the system. It travels at constant velocity 
vw and, as seen in Fig. [2j has a pronounced assymetric 
shape with a steep front and an extended tail. The an- 
harmonicity of this shape, together with previous results 
[26| , provide a simple explanation for the disconinuity of 
the global order parameter: The wave is born as a result 
of a linear instability and inherent noise. Once it exists, it 
grows to a large final size because nonlinear attenuation is 
ineffective. The definition of fl involves a spatial average 
which is dominated by the extended spatial region behind 
the peak of the wave, where local order is much stronger 
than in the corresponding homogeneous ordered state. 
Even though the area ahead of the wave front is mostly 
disordered with |w| ps 0, it cannot compensate the huge 
contributions to £1 from the densest part of the wave. As 
a result of this biased average, the global order parame- 
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FIG. 3. Order parameter versus noise r\ for several system 
sizes L. Parameters: R — 0.1, vo = 1. 

ter is much larger than in the homogeneous ordered state. 
These waves have remarkable properties. For example, 
Fig. QJb) shows that their speed vw is just slightly be- 
low the maximum possible speed Do, but larger than the 
speed of sound in the disordered phase, vs — v /V2- 
The waves are always supersonic with Mach numbers be- 
tween 1.19 and 1.41. The particles at the highest (that 
is densest) point of the wave are so strongly aligned that 
their average speed u = |w|/p is only slightly less than 
vw- Thus, most particles in the wave crest travel with 
the wave and do not just undergo restricted local motion. 
The particles just ahead of the wave front have low den- 
sity and display disordered motion. They do not "feel" 
the wave coming since it arrives faster than the speed of 
sound. Their territory is invaded and a fraction of them 
becomes strongly aligned and joins the wave for a while. 
This motivates the term invasion wave. 

In agreement with direct simulations 37|, I observe 
that the invasion wave has a perfectly straight front, 
perpendicular to its direction of motion. To acceler- 
ate the numerics, I take advantage of this apparent 
one-dimensional nature and drastically reduce the y- 
extension of the box to L y — 2, creating a very elon- 
gated simulation box. Examining large box lengths L x , 
one realizes that the maximum density in a wave and the 
steepness of the leading edge depends on system size in a 
very sensible way that transcends the traditional mean- 
ing of "finite size effects". In particular, the maximum 
density in the wave is proportional to L 2 , , and the width 
of the peak scales as 1/X|. In fact, the invasion wave can- 
not be seen as a localized perturbation of some mainly 
undisturbed medium. It is rather a global excitation of 
the entire system where, facilitated by periodic bound- 
ary conditions, all parts of the system are involved and 
particles everywhere adjust accordingly. In Fig [2th), by 
scaling the x-coordinate by L~ 2 and the density by L 2 X it 
is demonstrated that here is a master curve for the shape 
of an invasion wave. As a consequence, the maximum 
density gradient at the leading edge is proportional to 
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[3jthe hysteresis region grows with system size because the 
properties of the underlying soliton-like wave strongly de- 
pends on L. Compared to equilibrium systems, the phase 
behavior depicted in Fig. [3] looks unusual. Nevertheless, 
it semi- quan titatively agrees with direct simulations of 
theVM MM- 
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FIG. 4. Snapshots for the head-on collision of two soliton- 
like waves. Fhe sequence starts with two well seperated peaks 
close to the x-axis running towards collision with their steep 
fronts facing each other. At the latest time, the peaks are 
separated again after a successful "tunneling" through each 
other and now run towards the edges of the box. 



L\. To understand why these waves have not been de- 
tected by means of hydrodynamic approaches 25j, |30j] an 
effective Knudsen number, Kn = X/r is defined. Here, 
A = vqt is the mean free path of a particle and r is the 
radius of curvature of the density profile at the tip of the 
wave. It turns out that Kn is never smaller than about 
1/2, not even in the smallest systems that allow wave 
formation. At such Knudsen numbers, a structure with 
a small internal scale moves quite a large distance in one 
time step, causing hydrodynamic gradient expansions to 
diverge. 

Assuming a homogeneous ordered state, the mean-field 
theory of the VM 26j predicts that the flocking transi- 
tion is continuous. To check whether this is still true for 
inhomogeneous ground states, I calculate the order pa- 
rameter as a function of noise for different system sizes. 
According to Fig. [3J at the smallest size L — 4 where 
the system is still homogeneous, a continuous transition 
occurs. However, in bigger systems, very close to the 
predicted threshold r\c — 0.4361, ft still has a large value 
because of persistent waves. To accurately establish the 
order of this transition I used the following simulation 
protocol: An ordered state with a stable stationary wave 
at very low noise 770 is created. This inhomogeneous state 
serves as initial condition for a run with slightly larger 
noise n\ > r/Q. After convergence, the noise is increased 
again. This way, a sequence of stable inhomogeneous 
states with large order parameters is obtained, even at 
noise values a few percent above r\c- Approaching the 
threshold from the higher noise side confirms that the dis- 
ordered state with f2 = is stable for n > rjc- Hence, the 
flocking transition shows hysteresis; the wave state can 
coexist with the homogeneous state over a finite noise 
range An s=s 0.045?7c for L = 300. This concludes the 
proof for the discontinuity of the order-disorder transi- 
tion for large L, at the mean-field level. As seen in Fig. 



One of the defining properties of a soliton is the abil- 
ity to pass through each other without destruction. To 
perform this "soliton test" , I prepared stationary waves 
in two different systems with slightly different sizes, 
L x = 299 and L x = 300 and ensured they run in opposite 
directions. After the waves became stationary, I "glued" 
the two boxes together leading to a longer system with 
L x = 599. A series of snapshots of the time evolution of 
this two peak system is shown in Fig. @] At the earli- 
est time one sees two peaks running towards each other. 
Eventually, they start to overlap and form a large single 
peak. A while later, the two peaks reemerge with almost 
undisturbed shape like a conventional soliton. Watch- 
ing the time evolution through repeated collisions reveals 
that if the waves have a tiny height difference initially, 
this difference is amplified in every encounter. The bigger 
soliton takes a few particles away from the smaller one 
in every meeting until only one peak survives. Gradual 
coarsening also occurs inbetween collisions or when waves 
travel in the same direction. Therefore, true stationary 
states have only one peak. 



In conclusion, using kinetic theory I analyzed super- 
sonic waves that occur in the Vicsek model once the sys- 
tem exceeds a critical size. I demonstrated that the waves 
show hysteresis and alter the order of the flocking transi- 
tion from continuous to discontinuous. A phase diagram 
with atypical finite size effects was calculated. I argue 
that these waves were probably not detected in previ- 
ous analytical approaches because of their large Knud- 
sen number. This provides an explicit example of an 
active particle model where hydrodynamic equations fail 
to correctly describe one of the most important prop- 
erties - the nature of the phase transition. This could 
have implications for other, more sophisticated, models 
of active matter. While my calculations neglect correla- 
tions which are relevant at small particle speeds, they do 
demonstrate a novel mechanism to induce a first-order 
flocking transition. How relevant this mechanism is in 
the low speed regime and for particles with nonzero size 
remains an open question. I speculate that some aspects 
of this approach remain useful. For example, it might 
lay the ground for a theory in terms of interacting quasi- 
particles that represent soliton waves. 
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